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Abstract 

We investigate the particle and kinetic energy densities of harmonically trapped fermion gases at zero 
temperature in arbitrary dimensions. We derive analytically a differential equation connecting these 
densities, which so far have been proven only in one or two dimensions, and give other interesting 
relations involving several densities or the particle density alone. We show that in the asymptotic 
limit of large particle numbers, the densities go over into the semi-classical Thomas- Fermi (TF) 
densities. Hereby the Fermi energy to be used in the TF densities is identified uniquely. We derive 
an analytical expansion for the remaining oscillating parts and obtain very simple closed forms 
for the leading-order oscillating densities. Finally, we show that the simple TF functional relation 
Ttf[p] between kinetic and particle density is fulfilled also for the asymptotic quantum densities 
r(r) and p(r) including their leading-order oscillating terms. 

1 Introduction 

The recent experimental success in confining degenerate Fermi gases in magnetic traps [GJ has 
triggered a renewed interest in the theoretical description of the particle and kinetic energy densities 
of harmonically trapped fermions, both at zero @, |, |, |, |, 0, | and finite temperatures g, 



10[. Whereas most of these investigations were limited to d = 1, 2 or 3 dimensional systems, 
general analytical expressions for the particle and kinetic energy densities valid for any dimension 
d were given in 0, Il0|. The present paper is devoted to a discussion of exact and asymptotic 
results for particle and kinetic energy densities in arbitrary dimensions. We put an emphasis on 
separating their smooth and oscillating parts, identifying analytically the smooth parts with those 
obtained in the semi-classical Thomas- Fermi (TF) limit, and finding simple closed expressions for 
the asymptotically leading oscillating parts. Using the latter, we find that the TF functional relation 
t tf[p] between kinetic and particle densities holds also for the densities including the leading-order 
quantum oscillations, confirming a recent numerical observation Q of this surprising result. 

We start by recapitulating some basic definitions given in We consider a system of N 
spin-half fermions, trapped in a spherical harmonic potential in d dimensions 



2 

V(r) = ^-r 2 , r 2 = x( + x z 2 + ... + x z d . (1) 



We assume that the lowest M + 1 shells are filled, so that the ground state is non-degenerate. The 
local density then depends only onthe radial variable r and is given in terms of the wavefunctions 
&0O by 

p(r) = 2£ |^(r)| 2 , (2) 

€i<E F 
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where the Fermi energy Ep can be identified with the highest filled level, Ep = cm = M + d/2, and 
a spin degeneracy factor of 2 has been included. For the kinetic energy density, we consider the two 
equivalent expressions 

r(r) = ~^-2j2 <A*(r)V 2 <Mr), (3) 



"if/") |^2^ |V^(r)| 3 . (4) 

[Note that in the standard literature on density functional theory, r(r) sometimes denotes the 
quantity which we here call Ti(r).] In the presence of time-reversal symmetry (which implies the 
particle number TV" to be even), they are simply related by 

1 h 2 

r(r)=n(r)- - — V 2 p(r). (5) 



A convenient quantity is their average 



e(r) = i[r(r)+n(r)] > (6) 



which is numerically known M, 11] to be a smooth function whereas r(r) and T\(r) have oscillations 



that are opposite in phase (see also figure ^ below). We can then express r(r) and T\{r) in terms of 
£(r) and V 2 p(r): 

rir) = Or) V 2 p(r) , n(r) = £(r) + i ^ V 2 p(r) . (7) 

Eqs. (^) - ([/]) are exact expressions; for a given number of particles or filled shells in any dimension 
they can be computed exactly using the known quantum-mechanical wave functions. 

It is well-known numerically that for large numbers of particles, the exact densities go over to 
the smooth TF densities, as also demonstrated in the two following figures. In figure [j] we compare 
exact particle densities p(r) with those of the TF model (see section |3] for their precise definition) 
for two particular cases in two and three dimensions. In figure we show a comparison between the 
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Figure 1: Particle densities of fermions in harmonic traps. Left panel: N = 380 particles filling 
19 shells (M = 18) in 2D; right panel: N = 440 particles filling 10 shells (M = 9) in 3D. Solid 
lines: exact quantum-mechanical densities p(r), dashed lines: Thomas-Fermi densities pri?(r). 
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Figure 2: Kinetic energy densities of N = 110 fermions (M = 9) in a 2D harmonic trap. Left 
panel: three equivalent quantum-mechanical densities integrating to the exact kinetic energy. 
Right panel: £(r) and Thomas-Fermi density Trpir)- 

three forms @, @) and (||) of the kinetic energy density, along with that of the TF model. In all 
cases, the exact quantum-mechanical densities oscillate around the TF densities except near the 
turning points where the latter go to zero by definition. In particular, the two densities r(r) and 
7*1 (r) have oscillations with exactly opposite phase, which renders their average £(r) smooth and 
following very closely the TF kinetic energy density. 

This paper is devoted to examining various issues related to the correspondence between the 
exact quantum- mechanical densities and their TF limits, which in the framework of density func- 
tional theory correspond to the local density approximation (LDA). We first derive a homogeneous 
differential equation for the particle density p(r), linking it to the kinetic energy density £(r) in 
arbitrary dimensions. This has been one of the challenges in the density functional theory and has 
only been proved recently || in d = 2. An interesting consequence of this equation is that the Fermi 
energy to be used in the TF model is given analytically in terms of the number M of the highest 
filled shell. 

The surprising success of the TF or LDA formalism in reproducing the averaged quantum- 
mechanical results nearly all the way down to the turning point raises an important question about 
the domain of validity of the TF approximation. We examine this question in detail. Using a Taylor 
expansion of the exact densities we show that the leading-order terms in the limit of large particle 
numbers or, equivalently, large Fermi energies reproduce exactly the TF densities. Furthermore, 
both smooth and oscillating corrections to the TF limits can be gleaned from the Taylor expansions. 
We also use a calculation of the moments of the densities to buttress our conclusions. While the 
smooth parts may easily be given exactly for even dimensions, they are harder to obtain in odd 
dimensions, necessitating the asymptotic Stirling expansion of Gamma functions. Nevertheless, 
the leading-order results and their corrections provide clear pointers to an understanding of the 
numerical results. In particular, we succeed in extracting the asymptotically leading-order oscillating 
terms of the densities in a very simple closed analytical form valid for any dimension. The limitations 
of these asymptotic terms, and interesting relations between them, are investigated numerically. 

The remainder part of this paper is organized as follows. In section ^| we give explicit expressions 
for the densities using the definitions given above. Some of the results are already contained in g 
but are simplified here and given explicitly, in order to make the further analytical calculations more 
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transparent. Section |3| is devoted to the derivation of exact differential equations, valid in arbitrary 
dimensions, that connect the kinetic energy and particle densities, and of an equivalent integro- 
differential equation for the particle density that has the form of a Schrodinger-type eigenvalue 
equation. In section [|] we discuss the asymptotic limits of the exact densities for large Fermi 
energies A and discuss their precise relation to the TF densities. What emerges clearly from this 
analysis is that the densities can be uniquely separated into smooth and oscillating parts. The two 
leading-order terms in A of the smooth parts constitute analytically the TF densities, whereas the 
oscillating parts are of relative order 1/A (or lower), so that they vanish in the asymptotic limit. In 
section [| we discuss the oscillating parts of the exact densities. In particular, we show that their 
leading-order contributions, which for any dimension can be given analytically in terms of Bessel 
functions, are eigenfunctions of the radial Laplace operator. The TF functional ttf[p\ is shown to 
hold for the asymptotic densities including the leading-order oscillating terms. The full oscillating 
parts of the densities can be expanded as series of higher-order Bessel functions, which we exhibit 
analytically for even dimensions. Some numerical tests of our asymptotic relations are presented in 
section bA. In section || we finally present a summary of all the new results obtained in this paper. 
Some technical details are given in three appendices. 



2 Analytical densities for filled spherical shells 

Following the general definitions given in the previous section, we quote here the exact expressions 
for the densities derived in and give some new and more general analytical expressions which 
serve as starting points for the following investigations. In the remainder of this paper, we work 
throughout with dimensionless units corresponding to the choice m = u = K = 1. The exact 
quantum-mechanical densities p(r) and £(r) can be written Q, for a given number M of the last 
filled shell, in the simple form 

, M 

p{r,M) = _2^4? ) _ At (-irL M (2x)e^, (8) 

[i=0 
1 d M 

£(r,M) = —-Y^G { S_^-lYL,{2x)e-\ (9) 

/i=0 

where x = r 2 and L ti (2x) are Laguerre polynomials. The coefficients F„ , G„ are given by 

[u/2] [„/2] 

= v + 1 + iy + 1 - 2m) gW , = (y + l) 2 + £ (" + 1 " 2™) 2 9$ , (10) 

m=l m=l 

where [u/2] = integer (u/2). The factors gffl are defined through 

fl _ x) (l-d/2) = 1 , (d) m (d) = J_ r(d/2 + m-l) 

[ X) 9m ml r(d/2-l) ' 1 j 

(2) 

The results are particularly simple for d = 2, where all the gin are identically zero: 

FP = v + 1 , G( 2 ) = {u+ l) 2 . (12) 

In Q it was stated that F^ and Gv can be given in closed analytical form only for even 
dimensions. We now have found [12] that this actually can be done in any dimension d, albeit 
separately for even and odd indices: 

p( d ) r , /2 | nd/2 + n) {d) r(d/2 + n + l) 

F 2n - (d/2 + 2n) n , r(d/2 + 1) , ^ 2 „ + i - 2 „ !r(d/2 + 1) . (13) 
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Prom (|l^) one finds that the obey the following summation relations: 

M M 
2 ^- 1 )M-^ F (d) =F W +i (Meven), 2 £(-l) M ~^W = Fj? (M odd) . (15) 

Using the above formulae, we can readily give closed expressions for the densities at r = 0. For 
even and odd M, respectively, we get for the particle densities 

, 2 r(M/2 + d/2 + l) 2 (M + d)\\ 

/Wo,m) = ^ r(d/2 + 1)r(M/2 + 1) =^^?r (Meven) ' 

2 r(M/2 + d/2 + 1/2) 2 (M-l + d)!! . 

Podd(0, M) = r(d/2 + 1)r(Af/2 + 1/2) = ^72 d !!(Af-l)!! (M ° dd) ' (16) 

and for the kinetic energy densities 

U*(0,M) = (4 ^// 2 | 2) Pevcn(0,M), 

Ud(0,M) = (4M ( + + 3 ^ +2) Pod d(0,M). (17) 

On the right-hand sides in ([!(]), which are easier to evaluate manually, we have used the double 
factorials defined (for integer n) by 

(2n)H = 2n(2n - 2) • • • 2 = 0Fr(2n + l)/2 n r(n + 1/2) , 0!! = 1, 
(2n-l)!! = (2n-l)(2n-3)---l = 2 n r(n + l/2)/0F, (-1)!! = 1. (18) 

Note that the densities (|0]) for any even M and the next odd M+ 1 are identical. This is due to the 
alternating parities of the quantum states in successive shells. While filling a shell with odd parity 
(ie odd M + 1), no contribution is added to the density p even (r, M) at r = 0. For even dimensions 
d, the expressions (ll~6|) and (|l~7| ) reduce to polynomials of order d/2 and d/2 + 1 in M, respectively. 
For odd d, they can be simplified only in the large-M limit discussed in section ||. 

For later reference, we give here also the total particle number N and energy E as functions of 
the shell number M, making use of the fact that the degeneracy of the m-th level e m = m + d/2 is 
given by the binomial coefficient ( m l_]~ 1 ) : 

E(M) = 2 J] ( m + /_ - ^ (m + d/2) = d (2M + d + D ^ + ■ (20) 
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3 Exact differential equations for densities in arbitrary 
dimensions 

As stated by Minguzzi et al H, one of the challenges in density functional theory is to be able to 
directly calculate the particle density, given the potential energy, without recourse to solving the 
Schrodinger equation for the wave functions. They have achieved this aim for harmonic confinement 
of independent fermions in d = 2. The central result of their paper Q is the following differential 
equation for the density 



e(r) 



i Ap(r) + p(r) (m + | - l - r 



(21) 



where A is the Laplacian operator A = V 2 in two dimensions. 

Indeed, we now prove that such a differential equation may be derived in arbitrary dimension 
d. Furthermore, using a "differential virial theorem" derived recently [^] for any d, we show that 
a linear, homogeneous third-order differential equation for the particle density may be established 
in arbitrary dimensions, like it has been done in the case of d = 2. 

In order to proceed, we note that the radial part of the Laplacian in d dimensions is given by: 



A 



d 2 (d - 1) d 
dr 2 r dr 



For the following derivation it is useful to express it in terms of the variable x = r 2 : 

. d 2 , d 

A = 4x—, + 2d — . 
ax z ax 

The Laplacian acting on the particle density (||) is given by 

M 



Ap(x) 



1 



2 J2 F$-v(-m(±c -2d- 8p)L^2x) + (Ad - 8)L'^2x)} e 



7^/2 



A»=0 



where we have made use of the differential equation for the Laguerre polynomials 

yL»(y) + (l-»)L;(y)+^(y)=0 ) 



(22) 



(23) 



(24) 



(25) 



and the primes denote derivatives with respect to the full argument of a function. Manipulating 
the sums in (\24\i, we obtain 



1 4 

- Ap(x) + (M + 1 + d/4 — x/2) p(x) = - £(x) + H{x) 
o a 



(26) 



where 1Z(x) is given by 
M 



= E(-!) M {[(M + l-p) - Gg_J L,{2x) + (d/2 - 1) F^L'^x)} e~\ (27) 



Ai=0 



Note that 1Z(x) is identically zero for d = 2, due to (|l^), so that (|2lf ) reduces readily to (^|). For 
arbitrary dimensions d, now, 72. (x) can be cast into the following form (see appendix [A| for details): 



n(x) 



p{x) 



d " v " / 4 

Substituting this into (p6[), we arrive at the differential equation 



(d + 2) 



1 



Ap(r) +p(r) A 



1 



(28) 



(29) 
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which is valid for any d and obviously reduces to ( pl| ) for d = 2 as given in || . Here Am is given by 



\ M = M+(d+l)/2. 



(30) 



which corresponds to the average of the highest filled and the lowest unfilled level and will be 
identified as the Fermi energy of the smooth densities in the large-M limit (see section [| below). 

Having established equation (29), we now proceed to derive equations for the particle density 
alone by eliminating the kinetic energy density £(r). To this end we use an integral equation which 
has recently been derived in [1C] for arbitrary dimensions d, and its equivalent differential form: 

d 



d f°° 
£(r) = — / r'p(r')dr' 
2 Jr 



~rp(r) 



(31) 



The differential form had earlier been derived for d = 2 in p| and termed "differential virial the- 
orem". Its version for d = 3 had been correctly guessed in 0. Inserting the integral form of £(r) 
into <M), we obtain the following integro-differential equation for p(r) with eigenvalue Am 



1 a / s 1 2 / x (d + 2) 
_ Ap(r) + - r 2 p(r) + 



r p(r) dr = Am p(r) 



(32) 



This has the form of a Schrodinger equation for p(r) with an additional non-local potential. Differ- 
entiating both sides of (32), we can write it as a third-order differential equation for p(r): 



1 d 



8 dr 



o ~T~ Ap(r) + [ Xm - ~ r — p(r) + -r p(r 



d 



dr 



0. 



(33) 



This equation has been previously derived for d = 1 in 0] and for d = 2 in ||. Its form for d = 3 
has been surmised and numerically tested in ||, and general solutions for p(r) in the case d = 3 
were discussed in Q. 

Using the relation (0), we may eliminate the Laplacian term in (|29| ) in favour of the kinetic 
energy density r(r) to obtain the following exact relation between three densities 



r(r) = p(r) (\ M ~\ r2 ) ~\ ^ 



(34) 



which will prove useful in the discussion of section [|. 

Integrating the differential equation (29) over the whole space, using the measure 



d d r(...) 



IT 



d/2 



T(d/2) Jo 



x d / 2 - x Ax (...), 



and exploiting the virial theorem, we find 

XmN(M) 



d+1 
d 



E(M) . 



(35) 



(36) 



which using (n% and (po|) unambiguously confirms the choice of the Fermi energy (pO). 
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4 Large-M limits and relation to the Thomas- Fermi 
densities 

We next address the question as to how the well-known Thomas-Fermi (TF) results are obtained 
from the exact densities (pi), (P) in the limit of large particle numbers or, equivalently, in the large-M 



limit. The TF densities are given by (see eg chapter 4 in [16]) 

PTF(r,X) = j^- dY ^(\-r^e(X-r"/2), (37) 

Mr,X) = r TF (r,X) = Jp~yp, (dT2) r^)^ ~~ r 2 /2) d l 2+1 6(A - r 2 /2) , (38) 

where A is a Fermi energy that is usually determined by normalizing pTF( r , A) to the exact particle 
number N. The Heavyside step function 0(A — r 2 /2) ensures that these densities are identically 
zero outside the classical turning point r$ = \/2A. 

The proof that the above TF densities can be obtained from the exact quantum-mechanical 
densities in the large-M limit is by no means trivial, since the summations in (||) and @ cannot 
be done analytically in closed form. To our knowledge, a rigorous proof for d = 2 has been given 
only very recently in JL0|J , whereby A was found to be A = M + 3/2. In Q it was shown for d = 1 
that the exact density at r = is asymptotically given by ptf(O) for large particle numbers, and 
an expression for the leading-order oscillating part Sp(x) = p(x) — ptf(x), valid for small enough 
x, was given. 



We first observe that the differential equation (2£) is also satisfied by the corresponding TF 
densities in the limit of large M, and hence of large \m, since the Laplacian acting on the TF 
density produces terms of lower order in Xm- In the following we shall show, indeed, that for 
any dimension d the exact densities given in (^|) and (^) in the large-M limit go over into the TF 
functions (|37|) and (|38[), plus smooth corrections of order 1/M 2 relative to the leading-order TF 
terms, plus oscillating terms that are of relative order 1/M or lower and will be discussed in section 
[|. Hereby it is essential that the Fermi energy A be identified with the Am given in (|30|). Whereas 
the proof given in |l0[ for d = 2 can easily be generalized to any even d (see section |5] and appendix 
|C|), the cases with odd d turn out to be more difficult to handle. We will first discuss the separation 
of smooth and oscillating parts in general. In the two ensuing subsections, we shall then use Taylor 
expansions of the exact densities and a calculation of their moments in order to prove our claims. 

4.1 Separation of smooth and oscillating terms 

For the following discussions, it is useful to separate the densities into average and oscillating parts: 

p(r) = p(r) + 5p(r) , £(r) = £(r) + 5£(r) . (39) 
This can be achieved using the following trick. We define the smooth parts by 

p{r,M) = -[p even (r,M)+p odd (r,M)], 

f(r,M) = i[£even(r,M)+£odd(r,M)], (40) 

whereby the formal expressions obtained from (||) and @ for even and odd M must be taken at the 
same fixed value of M and inserted on the right-hand sides above, irrespectively of M being odd or 
even. In the same fashion, the oscillating parts are defined by 

5p(r,M) = (-l) M ±[p eyeil (r,M)-p odd (r,M)], 

6£(r,M) = (-l) M ^[£even(r,M)-£odd(r,M)]. (41) 
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The resulting expressions are then valid for either even or odd M. The sign (— 1) M in front of the 
oscillating terms is easy to understand due to the alternating parities of successive oscillator shells, 
as discussed in section Q after the equations (|l^) for the density at r = 0. The alternating sign of 
5p(Q, M) for successive shell numbers M is intimately connected with the oscillating behaviour of 
5p(r, M) as a function of r. For even M, it must have a maximum at r = since an even shell always 
contains an s (zero angular momentum) state. After filling the next (odd) shell, it has acquired a 
minimum, and so on. This is also illustrated in the numerical examples shown in figure |l[ 

From the explicit results given below it can be verified that inserting (^) into (|29|), the smooth 
parts cancel identically and the following differential equation holds for the oscillating parts alone: 

- A5p(r) + 5p(r) (^Xm — — r 2 ^j 



(d + 2 

Similarly, all smooth terms cancel also from (|34|), so that we have the equation: 



(42) 



5r(r) = 6p(r) (\ M - i r 2 ) - 1 <5£(r) . (43) 

We finally note that the differential equation given in (^|) for the full densities holds separately for 
their smooth and oscillating parts as well as for the TF densities alone: 

~j~£( r ) = ~^ r P( r ) ' dr^^ = ~^ r5p ^ ' ^^,TF(r) = p TF (r) . (44) 

Let us now look at the densities at the centre of the system, r = 0. As stated at the end of 
section |2[ for even dimensions d the expressions (16) and ( |T7| ) for p(0) and £(0) are polynomials in 



M of order d/2. It is easy to see that the coefficients of their highest power in M are identical with 
those of Ptf(0, Xm) and £tf(0, Am), respectively, so that the latter become exact in the large-M 
limit. For odd d, we have to use the Stirling expansion of the Gamma functions in fll6|), valid for 
large M (cf p^| , eq 6.1.37), in order to find the same result. It is even more instructive to look 
separately at the smooth and oscillating parts. We find that in the large-M limit, the smooth parts 
reproduce the two highest powers of M in the corresponding TF densities at r = correctly: 

"<°' m » - (2^4 <2[1+0(l/M2>+ - 
R**> - (^(^2)4 Ari[1+0(1/M2)+ '-' ] - <45> 

whereby the series in 1/M stop at order M° for even d. The oscillating parts are found to be 

Sp(0,M) = H) M ^4 A!! " |1 + 0(1/M) + '-' 1 ' 

5mM) = H) M ^||Ar[l + (l/M) + ...]. (46) 

We see that 5p(0) is of order 1/Xm lower than the leading power in the smooth density, whereas 
the suppression factor in the <5£(0) is l/X^. This hints at the numerically known fact that £(r) is 
much smoother than p{r). 

4.2 Taylor expansion of the densities 

We now want to study how the above results extend to x > 0. In the appendix [B|, we give the 
explicit Taylor expansions of the exact densities (||) and (^) in powers of x = r 2 : 

oo oo 

p^ (x, M) = J2 P {m) (d, M) x m , {x, M)=J2 £ (m) (d, M) x m . (47) 

m=0 m=0 
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The expansion coefficients cannot be given in closed form for any d, but must be evaluated separately 
for even and odd values of M. Prom these, one can reconstruct separately the expansion coefficients 
of the smooth and oscillating parts according to (40) and (pl|). 

As shown in detail in appendix |b|, we find the following structure of the results for even dimen- 
sions d. The smooth parts are finite polynomials in x which have exactly the structure of the TF 
densities ( |37[ ) and ( |38| ) , including the step functions that cut these densities at the classical turning 
point. However, only the two highest powers of M appearing in the coefficients of each power of 
x agree with those in the corresponding powers of the Fermi energy Am = M + (d + l)/2 ( |30D 
appearing in the TF densities. In other words, the smooth parts go in the large-M limit over into 
the TF densities plus smooth terms of order 1/M 2 and lower relative to the leading TF terms: 



p{r,M) = P TF(r,X M ) [l + 0{M~ 2 ) + ... 
£(r,M) = StfMm) \l + 0(M- 2 ) + ... 



(48) 



which is the straightforward extension of (45). To put this structure into evidence, we list here 
the complete smooth densities for the lowest even dimensions, re-expressed in terms of the Fermi 
energy. For d = 2 we get, with A = M + 3/2: 



p )(r) = 
For d = 4 we get, with A = M + 5/2: 



— (A — r 2 /2) 0(A - r 2 /2) 

1 

2^ 



(A - r 2 /2) 2 - 1/4 G(A - r 2 /2) . 



(49) 



£( 4 )(r) = (A-r72) a -7A/4 + 9r 2 /8 9(A - r 2 /2) . 



r) = ^ [(A - r 2 /2) 2 - 3/4] 9(A - r 2 /2) , 



6ir 2 



(50) 



The oscillating parts, characterized by the alternating sign (— 1) , contain infinite power series 
in x that can be resummed into oscillating functions of x, as shown in appendix [C] and further 
discussed in section [|. After resummation, they are found to be of order 1/M or lower, relative to 
the leading-order terms of the smooth parts, so that they become negligible in the large-M limit. 

The situation is more complicated for odd dimensions d. This becomes evident when we compare 
the powers of tt appearing in the denominators of the exact densities (||), @ and their TF expressions 
(|37|), (|38|): whereas they are the same for even d, they are different by a factor y/ir for odd d. The 
reason is that in order to obtain the TF densities in the large-M limit, one has to perform a large- 
M Stirling expansion of the T functions in (|16|), as seen above for p(0). In doing so for specific 
odd values of d, we find that the resulting densities have the same structure as for even d. Their 
smooth parts yield the Taylor expansions of the TF densities (|37|), which here are infinite due to 
the half-integer power d/2: 



p^(r,M) = p^(0,A) 



(d) 



d x d(d-2) 



2 2A 



+ 



x 

2A 



1 + C(M~ 



+ ... 



(51) 



with A = \m given in (|30|), and analogously for the average kinetic energy densities ^ d \r,M). The 
oscillating parts are again at least of order 1/M relative to the leading smooth parts and therefore 
vanish in the large-M limit. 

Summarizing up to this point, we find that the densities p(r) and £(r) can uniquely be decom- 
posed into smooth and oscillating parts. In the large-M limit, the smooth parts p(r) and £(r) reach 
their TF forms 



, ( p8|) in terms of the Fermi energy Am given in (3C), plus terms of relative order 
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1/M 2 and lower at each power of x in their Taylor expansions. The oscillating parts 5p(r) and 8£(r) 
are of order 1/M and 1/M 3 , respectively, compared to the leading-order TF terms, as discussed 
explicitly in section |5[ This clearly establishes the TF densities as the asymptotic large-M limits of 
the exact quantum-mechanical densities. 



4.3 Calculation of moments of the densities 

As an alternative way to demonstrate the asymptotic equality of the exact and TF densities, we will 
now calculate their moments and show that they are all identical to the two leading orders of M. 
Since the densities are functions only of x = r 2 , it is sufficient to compute only the even moments 
of r, i.e., the moments of the variable x. We define them as 

C {m) = J d d rx m p ( x j } D {rn) = J ^ ^.m ^ (52) 

Integrating Z)( m ) by parts using ([H]) and (|35"|), one easily sees that the following relation holds: 

D {m) = d (m+1) 

(4m + 2d) K J 

Note that this relation for m = yields the virial theorem, since is the kinetic energy and 
twice the potential energy. The moments of the TF density are readily found to be 

Mm) [ Ad m ( x x T{d/2 + m) 2 m+1 {m+d) 

C TF =Jdrx P tf(x, Am) = r((j/2) X m ■ ( 54 ) 

The corresponding moments of the exact density might be obtained by direct integration of x m over 
the expression which leads to 

Cff = r( r / ( ^ ) m) 2 J2F^(-irF(-p, m + d/2; 1; 2) , (55) 

where F(a, b; c; z) is a hypergeometric function. However, we were not able to do the above sums 
analytically for arbitrary d. We therefore turn to the case d = 1, where the moments are more 
easily obtained from a summation of the corresponding matrix elements. The matrix element of x m 
in the n-th quantum state of the one-dimensional harmonic oscillator is found to be 

(n\x m \n) = |^ (-l) n F(-n, m + 1; 1; 2) . (56) 



The exact moments are therefore 

" 2 2m ml 



c $ = iS 2 E ™ + 1; i; 2) . (57) 



n=0 

The above hypergeometric functions are easy to create recursively by the relation 

(m + 1) F m+1 (n) = (2n + 1) JP" m (n) + mF m ^i(n) , (58) 

where 

^ m (n) = (-l) n F(-n,m + l;l;2), (59) 
and using F(-n, 1; 1; 2) = (-1)" (see ||, eq 9.121.1). This yields 

^bH = i, 

^i(n) = 2n + l, 
jT 2 ( n ) = 2n 2 + 2n + l, 

F 3 (n) = (2n + l)(2n 2 + 2n + 3)/3, (60) 
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and so on, after which the summation in ( |57[ ) becomes elementary. Below we give in table [T] (left 
part) a list of the first few moments, obtained in the TF approximation from (|54| ) and exactly from 
(p7|) ; expressed in terms of the Fermi energy A = M + 1. We see that the TF moments reproduce 
the exact ones up to terms of relative order 1/A 2 and lower, which are equivalent to terms of order 
1/M 2 and lower. 

The same result can, with some more algebraic effort, be found for higher dimensions by manual 
evaluation of (p5|). For even d, the moments are more easily obtained from a Taylor expansion of 
the Laplace transforms P(s) and X(s) of the exact densities, given in appendix Q in powers of s. 
For d = 2, the simplicity of (|l2|) allowed us to derive the recurrence relation 

' m + V 



qm 



1) 



m + 3 



2 A Cjg$ + m 



2 (~t(m- 



1) 



(d 



(61) 



In the rightmost column of table [T] we give the lowest quantum- mechanical moments for d = 2, 
and in table we give them for d = 3 and d = 4. Their leading terms reproduce in all cases the 
TF values (|54|). Note that the polynomials seen in these tables, translated into polynomials in M, 
appear exactly in the coefficients of the Taylor expansions of the densities for even d. This can be 
checked for d = 2 and d = 4 with the expansions given in the appendix |b|. 

This establishes the desired result that in the asymptotic limit, all moments of the quantum- 
mechanical density p(r,M) are identical with those of the TF density Ptf{t,\m) up to relative 
terms of order 1/M 2 or 1/\ 2 M . By virtue of (|53"|), the same conclusion can be drawn for the kinetic 
energy density £(r) since the proportionality constant in this equation does not depend on M or A. 



4f ] id = 1) 



Cj$ (d = 1 , A = M + 1) 



(d = 2 , A = M + 3/2) 



in 




1 

2 
3 
4 
5 
6 



2A 
A 2 
A 3 

5 \4 
4 A 

T A 6 
f A 7 



2A 
A 2 

A 3 + ±A 
f (A 4 + 2A 2 ^ 



| (A 5 + 5A 3 + |A) 

21 ( A 6 + 1QA 4 + 23 A 2 } 

f ( A 7 + f A 5 + 49A 3 + 



A) 



y 

2 
3 
2 
3 

4 

5 V" 1 2' 

90 \ 

l(A 7 + fA 5 + f A 3 -%X) 

16 (A 8 + 2U 6 + 567 A 4 



(A 3 - |A) 
(A 4 + iA 2 - 
(A 5 + | A 3 - 

a ( a6 + f a4 



t|a) 
1 Ha 2 



+ T6 A 



1575 ' 
256 i 



Table 1: The lowest moments of the density of the 1-dimensional harmonic oscillator 

with M + 1 filled shells, expressed in terms of the Fermi energy A = M + 1. Column 2: 
quantum- mechanical values from (57); column 3: TF values from (|5"4|). Column 4 contains 
the quantum-mechanical moments for the 2-dimensional oscillator; the TF values (not shown 
explicitly) are given by the leading term for each m. 



m 


C}$ (d = 3 , A = M + 2) 


C ( q Z ] (d = 4, A = M + 5/2) 



1 

2 
3 
4 


\ (A 3 - A) 
1 (A 4 - A 2 ) 
1(A 5 -A) 

24 (A + 2 A - 2"A ) 
| (A 7 + 7A 5_ 7 A 3_9 A) 


1 f\4 5 \2 , 9 \ 

IB" (A - 2 A + TeA) 

l/\6 5 \4 41 \2 , 45 \ 
15^ ~ 4 A ~ T6 A +64^ 

8 ( \7 1 7 \5 161 \3 , 153 \\ 
T05 ^ + 4 A ~ T6" A + 

2 r\8 i 7\6 133 \ 4 177\2 , 945 \ 
m {A + (A - —A - -j^-A + 256 J 



Table 2: Same as in column 4 of table [l], but for d = 3 and d = A. 
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5 Oscillating parts of the densities 

In this section we focus on the oscillating parts 5p(r) and <5£(r) of the densities as defined in (|4l|). We 
first derive series expansions of these oscillating densities and then exhibit that their asymptotically 
leading terms are of order 1/M and 1/M 3 , respectively, relative to the leading-order terms of their 
smooth parts. 

5.1 Expansions of the exact oscillating parts 

We employ here a method proposed recently in In terms of the variable t = 2x = 2r 2 , we 

define the Laplace transforms 

P{s)= e- st p(t)dt, X(s)= e- st £(t)dt, (62) 
Jo Jo 

For even dimensions d these can be evaluated analytically, as shown in appendix |C], and the trans- 
forms of the smooth and oscillating parts can be readily identified. Taking the inverse Laplace 
transforms of the smooth parts yields exactly the results discussed in the previous section and sum- 
marized in (|48|) , including the step function which puts the TF densities identically to zero at and 
beyond the classical turning point. The inverse Laplace transforms of the oscillating parts yield, 
after a suitable expansion for large M which is shown in appendix]^, the following series for d = 2: 

f 1 A 1 1 

6p m {r) = L^^ Mz) __ x 2j 2{z) __ x sj 3{z) __ x 4 Mz)+ ^^ (63) 

(—-\\ M ( 1 A 1 1 

^ {2) W = ^—^! [ -2xJ l (z) + -x 3 J 3 (z) + -x 4 Mz) + m x 5 J 5 (z) + ..j, (64) 

with A = M + 3/2 and 

z = 2 V2\x~ = 2 rV2X . (65) 
Hereby we have defined the functions J^{z) by 

Mz) = T(p + 1) (£f (M > -1) (66) 

in terms of the standard Bessel functions Ju(z). Note that the Jp[z) functions are eigenfunctions 
of the radial Laplace operator, as discussed further below, and normalized such that ^(0) = 1. 
For rf = 4 we find the results 

*P (A) (r) = [ -^{\Mz) + -xMz) + -x 2 J 2 (z)--x*J 3 (z)--x 4 Mz) + ..j, (67) 

^ (4) W = [ -^{Mz)-- 4 x 2 J 2 (z)--x 3 J 3 (z) + -x 4 Mz) + m x 5 J 5 (z) + ..j, (6^ 

with A = M + 5/2. It is readily checked that in both cases, the above results satisfy the differential 
equation in ( |44| ) connecting 5£(r) with 5p(r). 

Using the known Taylor expansions of the Bessel functions, the above results yield exactly the 
expansions of the oscillating terms given in appendix [B| and discussed in the previous section. For 
the following it is important to note that the two leading powers in M of each Taylor expansion 
coefficient are correctly reproduced by the leading-order Bessel functions in the above series; the 
higher-order terms just correct the terms of order 1/M 2 and lower in each coefficient. 

We were not able to find the same expansions of the oscillating densities for odd d, since 
their Fourier transforms could not be obtained in a tractable analytical form. We have contented 
ourselves with the observation that the leading-order Bessel functions, given below for arbitrary d, 
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are recovered from the Taylor expansions discussed in the previous section, again up to terms of 
relative order 1/M 2 and lower in each Taylor coefficient. 

The convergence of the above series is tested numerically in section |5.4| below. They are asymp- 
totic in nature and converge only for sufficiently small values of x, ie, sufficiently far from the 
classical turning points. At x = 0, only their leading terms (ie, the first Bessel function in each 
series) contribute, for which we give some particularly nice closed forms in the following subsection. 



5.2 Asymptotically leading oscillating parts 

We now limit ourselves to the leading-order terms of the oscillating densities, which we will denote 
by the subscript "as" for "asymptotic" and which are given by the first Bessel functions in each of 
the above exact expansions. For arbitrary d and M they can be written in the following forms: 

WD - fe^(^)V), (69) 

with 



p = d/2-l, z = 2ry/2X^, X M = M + (d + 1)/2 . (71) 
Note that the above asymptotic forms fulfill the same differential equation as that given for the 



full oscillating parts in fl44|). Comparing to the TF densities in (p7|), (38), we see that Sp^^r) and 
<5^ as (r) are of relative order 1/Xm and 1/A|^, respectively. This confirms the numerical finding that 
even for moderate values of M, the densities £(r) appear to be smooth, while the oscillations in 
p(r) are clearly visible (see figures [l] and |2|). 

It is rather easy to verify that the asymptotic densities (]69|), ( |70|) are eigenfunctions of the 
Laplace operator for any p, and hence for any d, with eigenvalue — 8Xm- 

A 6 Pas (r) = -8X M 5 Pas (r) , A 5£ as (r) = -8A M ^ as (r) . (72) 



Note, however, that this is not true for the full oscillating densities given in fl63|) - fl68|) above. The 
equation for (5/? as (r) in ([72]) can also be obtained from the exact differential equation (^) derived 
in section ||. Taking the asymptotic limit of ([42;), retaining only the two leading powers in Xm, the 
5£(r) on the lhs and the potential term on the rhs can be neglected so that the first equation in 
7% immediately follows. 

Gleisberg et al || have studied the asymptotic expansion in the case d = 1 in which our result 



(—l) M 

6pW ( x ) = cos(2xv / 2A) , A = M + 1 , (73) 

7TV 2A 



for 5p as becomes 



where x here is the dimensionless coordinate. This result agrees with the oscillating part of their 
result (cf P|, eqs 16 and A9), if we replace our Fermi energy A by M + 1/2 which is that of the 
exact quantum system. (Note that these authors did not include the spin factor 2 in their density.) 

We emphasize that the asymptotic forms (|69|), ( f70| ) give the exact values of the oscillating 
densities at r = 0, since all higher correction terms in the series (|63|) - (|68|) are zero at the origin. 



As shown in the numerical investigations of section 5.4, they describe the quantum oscillations of 
the exact densities near the origin very precisely, where the densities are largest and perhaps most 
easily detected experimentally. 

Next we derive the leading oscillating part of the kinetic energy densities r(r) and t\(t) which 
are seen in figure § to exhibit much stronger oscillations than £(r). To derive the asymptotic forms 
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of their oscillating parts, we start from (fi3| ) and keep, as above, only the two leading-order terms 
in Xm on both sides of this equation. We then find the following simple relation 

tfT M (r) = \ M S Pas (r), (74) 

which will be tested numerically in section Using the relations (|5|) and (|72|), we finally obtain 

ST 1>as (r) = -X M 5pas(r) = -Sr^r) , (75) 

which confirms explicitly the numerically known fact (cf figure ||) that the oscillations in T\{r) are 
opposite in phase to those in r(r). 

5.3 Validity of the TF kinetic energy density functional 

In 0, the surprising observation was made that the TF functional for the kinetic energy density 



h 2 iird 
ttf[Ptf\ = - — 



2m (d + 2) 



d -r( d - 



4 V2 



2/d 

1+2 d 

Ptf > ( 76 ) 



is fulfilled to a high degree of precision also for the exact quantum-mechanical densities p(r) and 
r(r), except near the classical turning point and beyond it. Upon integrating, in the case d = 2 one 
even obtains the exact kinetic energy. 

Using the above asymptotic results it is easy to prove the validity of the TF functional ( f76| ) at 
the level of the leading-order asymptotic terms. We define the asymptotic total densities by 

Pas(r) = Ptf{t) + 5 p^ , T as (r) = T TF (r) + <5r as , (77) 

insert these into the TF relation fl76|), expand the power [p{r)\ 1+2 l d up to first order in Sp^, and 
neglect all terms of relative order 1/Xj^ and lower. We then find that the TF functional relation 
( |76D is, indeed, fulfilled also by the total asymptotic densities ( |77[ ) 

T TF [pUr)] = Tas(r) + 0(M- 2 ) , (78) 

which explains the results found numerically in [|]. 

The result (|7S] ) is highly nontrivial. In principle, the TF relation ttf[p\ holds only for systems 
with constant densities where the LDA is exact. For smoothly varying densities, gradient corrections 
are known to exist, which may be derived in the so-called extended TF model and have been 
successfully applied to various fermion systems (see eg chapter 4). Nevertheless, our result 
( |78D shows that the simple TF functional applies asymptotically for harmonically trapped fermions 
also when including the quantum oscillations in both p(r) and r(r). 

5.4 Numerical tests of the asymptotic relations 

We first test the asymptotic form (|6^) of the oscillating part of the density. In figure £3, it is compared 
to the full oscillating part 5p(r) defined by ( |4l| ) for iV = 420 particles in a 2D trap. The cusp in the 
exact 5p(r) reflects the cut-off of the smooth TF-like density p(r) in ([49|) at the classical turning 



point. We see that 5p as in (69) is correct only for small distances r from the centre of the system 
and completely fails near the classical turning point. The situation improves somewhat when we 
include higher terms of the expansion discussed above. In figure |] we show the same comparison 
for N = 1722 particles in 2D, with and without including the terms up to the order x^J%{z) in 



the Bessel function expansion (63). Whereas 5p as is only sufficient up to r ~ 1, the Bessel series 
truncated at order 8 becomes good up to r ~ 2.5; the phases of the oscillations are correct even 
further out. But clearly, an infinite number of terms would have to be included before reaching an 
agreement in the region of the classical turning point. 
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r 



Figure 3: Oscillating part of the density for N = 420 particles in a 2D harmonic trap (M = 20). 
The solid line gives the full oscillating part as defined by ( |4l| ) , and the dashed line its asymptotic 
form given by ( |69[ ) for \x = 0, A = 21.5. 




Figure 4: Same as figure |, but for N = 1722 particles (M = 40, A = 41.5). The dotted line 
gives the results obtained when including the terms up to O^J^z)] i n the expansion (|63|). 

Whereas this result is not very encouraging, showing that the asymptotic oscillating densities 
( |69|) are only useful far away from the classical turning points, a much larger region of validity is 
found for the relation ([74|) between the asymptotic densities Sr as (r) and 5p as (r). This is demon- 
strated in figure [|, where we have plotted the full quantities 5r{r) and \ 5p{r) (with A = 62) for 
the N = 79422 particles filling 61 shells in a 3D harmonic trap. The agreement is quite good until 
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10 



\8p(r) 




2 4 6 8 10 12 

r 

Figure 5: Oscillating parts of densities for N = 79422 particles in a 3D harmonic trap (M = 60, 
A = 62). The solid line gives Sr(r), and the dotted line the quantity X5p(r) according to the 
asymptotic relation (f74|) . 

rather far away from the centre; only near the turning point the proportionality breaks down. The 
same is shown in figure |6] for the N = 271502 particles filling 41 shells in a 4D harmonic trap (with 
A = 42.5). To set the scale, the density p{r) and its TF limit are also shown in the figure. The 
particle number is so high here that the oscillations in p(r) can only be recognized very close to the 
centre, which demonstrates once more the excellent approximation provided by the TF density. 
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Figure 6: Same as figure || but for N = 271502 particles in a 4D harmonic trap (M = 40) with 
A = 42.5. The dashed ancT dashed-dotted lines give the full density and its TF limit. 



Clearly, an appreciable amount of error cancellation must take place on the way to the relation 
(f74|) between the asymptotic oscillating densities. Its good validity is the basis of the performance 
of the TF kinetic energy functional ttf[p] including the quantum oscillations, as discussed above 
and established formally through (|7q). 
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6 Summary 



We have studied the particle densities and several equivalent forms of kinetic energy densities 
for harmonically trapped fermion gases in arbitrary dimensions d. We established various exact 
differential equations connecting these densities, and an equivalent integro-differential eigenvalue 
equation of Schrddinger type for the particle density. Some of these equations have been previously 
derived, but only in d = 1 or 2 dimensions. We have then studied in detail the asymptotic limits 
of the particle densities p(r) and the kinetic energy densities £(r) for large particle numbers and 
Fermi energies A. Using Taylor expansions, moments, and Laplace transforms of these densities, 
we were able to separate them uniquely into smooth and oscillating parts. The two leading-order 
powers in A of the smooth densities p(r), £(r) were shown analytically to go over into the well- 
known Thomas-Fermi (TF) densities, the remaining terms of order 1/A 2 and lower giving small 
corrections to the smooth parts. The Fermi energy hereby is given by A = M + (c£ + l)/2, where M 
is the quantum number of the highest occupied spherical shell. The oscillating parts Sp(r), 5£(r) 
are characterized by an oscillating overall sign (— 1) A/ that reflects the alternating parities of the 
successive shells. Their asymptotically leading terms could be given in simple closed expressions 
involving Bessel functions; they are eigenfunctions of the radial Laplace operator with eigenvalue 
— 8A. They are of order 1/A and 1/A 3 for 5p(r) and 6£(r), respectively, with respect to the leading 
terms in the smooth densities p(r) and £(r). The latter fact - the suppression factor 1/A 3 - explains 
the numerical experience that the kinetic energy densities £(r) are smooth even for moderate particle 
numbers. For even dimensions d, we have derived a formally exact expansion of the full oscillating 
densities in terms of Bessel functions, which was, however, found to converge rather slowly and only 
sufficiently far away from the classical turning points. Still, these results are useful for expressing 
the quantum oscillations near the centre of the system in closed analytical forms. Using our exact 
differential equations which hold separately for the smooth and oscillating parts of the densities, we 
could also establish the numerical fact that the kinetic energy densities r(r) and ri(r) are oscillating 
with opposite phases. An interesting proportionality between 5p(r) and 5r{r) was found to hold 
asymptotically everywhere except near the classical turning points. With our asymptotic relations, 
we could finally establish analytically the surprising recent observation Q| that the TF functional 
t tf[p] holds also when the quantum-mechanical oscillations are included in both r(r) and p(r). 

One of us (MB) acknowledges stimulating discussions with B van Zyl at early stages of these 
investigations. We are also grateful to R K Bhaduri for his continued interest. MB acknowledges 
the warm hospitality extended to him at the IMSC Chennai during a research visit, and financial 
support by the Deutsche Forschungsgemeinschaft. 
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A Evaluation or lZ(x) 

In this appendix we give details on the evaluation of the quantity IZ(x) defined in (^7j) and used in 
deriving the differential equation (29). Here we prove the following identity: 

n M 



K(x) 



H=0 

(d-2 



(d) 



G 



(<l) 



d 



i{x)-^—^- P {x). 



Proof: Using the identity 



n-l 



L' n (2x) = - ]T L m(1x) , 



(n > 0) 



(79) 



(80) 



m=0 



which may be derived using equations 8.971.2 and 8.974.3 of p5| , noting that L' (2x) = 0, and 
inverting the order of summation in the first line of (|79|), we get 



M 



fj,=0 



KM = ^7/2 E(" 1 ) M_M {\Q* + 1) F ^ ~ G^\L M ^(2x) 



(d-2] 



M-n-X 

■iff) £ L m (2x)J e ^. (81) 

m=0 



(The sum over m above must be put to zero when the upper limit of m takes the value —1). We 
may now reorder the double sum over the last term in such a way that the Laguerre polynomials 
can be taken out of the inner sum, leading to 



M 



K(x) 



7^/2 



J2(-l) M -' l LM-,(2x)H^e- x , 



(82) 



with the definition 



< } = (m + ±H d) - - (-1)" E(-ir^ 



(83) 



u=0 



For the sum over the Fu on the rhs above, we can use the formulae ( |15| ) given in section ^. 
Furthermore, from our explicit results (|l^) and ( |l~4|) one easily sees that 



(l + d/2)Gg +1 = (4n + 2 + d)F^ +1 , 



(l + d/2)Gg = (4r,+ l + d/2)F£) + (d/2-l)F£L 1 . 
Using these identities, it just takes some algebra to find the simple form 



(d-2) r 



G (d) _ F (d) 



(84) 



(85) 



valid for both even and odd fi. Substituting this into (|82| ) we finally obtain the desired result (79). 
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B Taylor expansions of exact densities 

We give here some Taylor expansions of the exact densities (||) and (||) in powers of x = r 2 . Their 
coefficients, defined in (f47|), cannot be given in closed form for any d and M, but only separately 
for even and odd M. The coefficients p^ (d, M) are simply the densities at r = given already in 
(16). For the next three coefficients of the density p(r), we find for even M: 



o( 2 ) 



(d, M) 



Pe7en(d, M) 



,(3) 



(2 + d + AM) 
2(d + 2) 



(86) 
(87) 



Pev'en(d,M) = 

and so on. For odd M we find: 

P ( l(d,M) = &(d,M), 



(d 2 + 6 d + 8 + 12 dM + 16M + 16M 2 



6(d + 4)(d + 2) 



P ( X(d,M), 



^odd 1 - 

Pol d (d,M) 



P { olM,M) 



(2 + 3d + 4Af) (o) 
2(d + 2) ^ddK^), 

(5 d 2 + 10 d + 8 + 20 dM + 16M + 16A/ 2 



&(d,M), 



(90) 
(91) 



6(d + 4)(d + 2) 

and so forth. Similar expressions are obtained for the coefficients ^ m \d, M). We can then recon- 
struct separately the expansion coefficients of the smooth and oscillating parts defined in (fiol) and 



(41). We give here the explicit expansions, valid for both even and odd M, for even dimensions, 
where the results look more transparent. For d = 2 we obtain 



1 

7T 



1 

2^ 

(-1) 



M + 



1 



(-1) 



M 



2tt 



1 - 2x \M + + x 2 (M 2 + 3M + 2 



(M 2 + 3M + 2 



x M + - + -x 



-x 3 { M 3 



1 



9 



13 



- AT + — M + 3 + 



(92) 



4tt 



Mr / 3\ 1 

x-x 2 (Af + -J+-x 3 ( ..1 A- + 3.1 / + 2 



2J 4 

2 



18 



M 3 + - M 2 + — M + 3 
2 2 



180 

For d = 4 we get explicitly 



+— x 5 (M 4 + 6M 3 + 14M 2 + 15M + 6) + 



(93) 



/4f( X ) - 715 



+ 



1 

-1) 

4^ 2 



M 2 + 5M + — 
2 



x I M + r ] + - x 



M 



M + 



x M + 5Af + 



11 



1 



15 



+ - x z AP + — AT + 17M + 



45 



+ 



1 

-1) 

8vr 2 



x 3 ( M 4 + 10M 3 + 35M 2 + 50M + 24) + . . . 



11 



M 



18 

M 3 + — M 2 + 17M + — 
2 4 



l-2x( M + - ) +x 2 f M 2 + 5M + 



(94) 



- x AT + 5Af H + - x z M + - 



+ 



-x 4 (M 4 
36 V 



11 



- x 



2 .3 ^ M 3 + ^ M 2 + 17M+ 45\ 



+ 10A/ 3 + 35AT + 50M + 24 + . . . 



4 J 
(95) 



Note that the smooth and oscillating parts in the above expansions of both densities are readily 
recognizable, the latter being multiplied by the overall sign (— 1) A/ . 
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C Asymptotic expansion of oscillating parts 



We start from the definition of the Laplace transforms (|62|) of the exact particle and kinetic energy 
densities (|8|) and (|9|). Using the integral (see Jl5| , equation 7.414.2) 

(s- 1/2)" 
(s + 1/2) 



/ e- st L n (t)e- t/2 dt 
Jo 



n+l 



(96) 



one obtains 



P (d) (s) 



M 



fi=0 



M-n 



1 d 



2 F\ d) 
=o 

M 



;i/2 



1/ 2 + s )m+i ' 

(1/2 - s y 



(97) 



We first discuss the simplest case, d = 2, and reproduce the results of |10| in a slightly modified 



form. Here, as for any even dimension d, the sums in ( |97| ) can be done analytically [12], with the 
results 
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8vrs 2 (1 + 1/2 S ; 
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(98) 
(99) 



These expressions fall into two parts which correspond exactly to the Laplace transforms of the 
smooth and oscillating parts of the densities defined in (|39|), whereby we recognize the alternating 
sign (— 1) M in front of the oscillating terms. The smooth parts give after Laplace inversion exactly 
the results (||). 

In the following we shall concentrate on the oscillating terms. We shall presently show that 
these terms yield oscillating functions of r (or x) and derive their systematic asymptotic expansion. 
To this purpose we extract a common factor appearing in both functions above 

E\ s ) - , - | ) (\ - ^ ) ■ A - M + .3/2 . ( 1 00) 



(1 + 1/2*)^+! U + l/2s 



Using the substitution s ■ 
expansion Jl^] for A — > oo: 



X/w we rewrite the first factor above and evaluate its asymptotic 
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(101) 

Next we re- insert w = X/s on the rhs above and perform a Taylor expansion of the second factor in 
( |100| ) in powers of l/2s. Note that this expansion converges since we can always deform the contour 
of the inverse Laplace transform, which is a vertical line to the right of the imaginary axis in the 
complex s plane, such that |l/2s| < 1. As a result, we obtain the following asymptotic expansion 
ofE(s): 

_1 X 1 X_ (>? 1_\ 1_ > 

8s 2 ~ 12s3 ~ 128s 4 ~ 480s 5 + 1 288 ~ 1024 J ~ 53670s 7 
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(102) 



Together with the prefactors appearing in (p8|) and (p9|), we can now Laplace invert the resulting 



series term by term. Hereby we use the general formula (see [13], equation 29.3.80) 



M2Vtx) 



(103) 
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where the Ju(z) are the cylindrical Bessel functions. Expressing the latter in terms of the 
defined in (]66|), we arrive at the results (|63|), fl64]) given in section [|. 

For the case d = 4 we just give here the Fourier transforms of the oscillating parts 



(1 - l/2s) M + 4 
(l + l^s) 1 ^ 1 



^) M _(-i) M (i-v^ 

d 8vr2 s (l + l/2 8 )^+i' ^04) 



which can be again be clearly separated from their smooth parts that yield the results ( |50| ) after 
Laplace inversion. Proceeding exactly as above, but now using A = M + 5/2, we arrive at the series 
(p7|) , ( |68|) of the oscillating densities given in section [|. 

Unfortunately, we could not find any tractable analytical expressions after the summations in 
(P7|) for odd dimensions d, which yield rather cumbersome combinations of hypergeometric series. 
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